今天跳tone一下,來說一下視覺化的東西,不然這幾天都貼了視覺化程式,卻一點說明都沒有。
基本上,一維的圖就不用特別了吧。重點是二維含地圖的視覺化。
主要可以用basemap跟cartopy兩個套件配合matplotlib完成二維含地圖底圖的視覺化。
這邊說明一下cartopy視覺化高解析可見光反照率資料-東亞
#讀檔的部分請參考第10天的內容,我這邊物件的名稱會維持一致
#得到投影資訊
print(sattif.GetProjection())
#得到起始角落點及結束角落點
#llx及lly分別為起始(通常是左上角)座標點,pixelx及pixely為沿著x及y軸之單位pixel距離。
#skx及sky為偏度,通常都是0。
#lrx,lry分別為結束(通常是右下角)座標點。
llx, pixelx, skx, lly, sky, pixely = sattif.GetGeoTransform()
lrx = llx + (sattif.RasterXSize * pixelx)
lry = lly + (sattif.RasterYSize* pixely)
上面程式會print出以下訊息
PROJCS["unnamed",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",128.5],PARAMETER["standard_parallel_1",30],PARAMETER["standard_parallel_2",60],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]
那可以看到這一串資訊內有一個PROJECTION["Lambert_Conformal_Conic_2SP"],代表這分資料的投影是使用蘭伯特投影,並且相關設定也一同提供。接下來就使用這些資訊畫圖
import matplotlib.pyplot as plt
import cartopy.crs as crs
import cartopy.io.shapereader as shpreader
from cartopy.feature import ShapelyFeature
import matplotlib.colors as mcolors
#先設定畫圖的底圖為蘭伯特投影地圖
#central_longitude對應的是central_meridian,central_latitude對應的是latitude_of_origin
#standard_parallels對應的非別為standard_parallel_1及standard_parallel_2
pj = crs.LambertConformal(central_longitude=128.5,central_latitude=0,
standard_parallels=(30,60))
fig , axs = plt.subplots(1,2,figsize=(8,8),subplot_kw={'projection': pj}) #設定底圖
axs[0].imshow(plt.imread("LCC_VIS_TRGB_2750-2022-09-24-10-10.jpg"))
axs[1].set_extent([llx,lrx,lry,lly],pj)
ishow = axs[1].imshow(albedo,cmap="Blues_r",zorder=1,extent=(llx,lrx,lry,lly))
axs[1].coastlines(zorder=2)
plt.colorbar(ishow,fraction=0.05)
plt.title(obsTime)
視覺化如下
可以看到圖雖然是正正方方的,但經緯度是斜的,這是因為現在視覺化的方式是將蘭伯特投影的資料放在蘭伯特地圖的底圖上,也就是同樣的角度看自己,當然會是正正方方的。
以下的程式,改成將蘭伯特投影的資料放在直角投影地圖的底圖上
fig1 , axs1 = plt.subplots(1,2,figsize=(8,8), \
subplot_kw={'projection': crs.PlateCarree()}) #設定底圖
axs1[0].imshow(plt.imread("LCC_VIS_TRGB_2750-2022-09-24-10-10.jpg"), \
extent=(llx,lrx,lry,lly),transform = pj)
axs1[0].set_extent([llx,lrx,lry,lly],pj)
axs1[1].set_extent([llx,lrx,lry,lly],pj)
ishow = axs1[1].imshow(albedo,cmap="Blues_r",extent=(llx,lrx,lry,lly),zorder=1,transform = pj)
axs1[1].coastlines(zorder=2 )
axs1[1].gridlines(draw_labels=True)
plt.colorbar(ishow,fraction=0.05)
plt.title(obsTime)
視覺化如下
可以看到經緯度變成直角了,但資料的呈現變成不是方正了,這是因為地球是橢圓的,所以這樣的呈現方式是維持橢圓資料的形狀放在以直角座標攤開的地圖上。
所以呢,只要掌握了地圖頭投影的正確概念,就可以畫出很多不同角度呈現資料的視覺畫圖型。